* hellwig-makeplots.do
* This file uses the simulated beta coefficients to make hellwig plots. Also does delta method
* -----------------------------------------------------------------------------
cd "~/Dropbox/Benton-Jordan-Philips/final-JOP-replication"

* some globals to simplify code below:
global myreshape "reshape long sigma2true_ sigma2_ll95_p sigma2_ll90_p sigma2_ll75_p sigma2_med_p sigma2_ul75_p sigma2_ul90_p sigma2_ul95_p sigma2_ll95_m sigma2_ll90_m sigma2_ll75_m sigma2_med_m sigma2_ul75_m sigma2_ul90_m sigma2_ul95_m sigma2_ll95_r sigma2_ll90_r sigma2_ll75_r sigma2_med_r sigma2_ul75_r sigma2_ul90_r sigma2_ul95_r sigma2_SDul_p sigma2_SDll_p sigma2_SDul_m sigma2_SDll_m sigma2_SDul_r sigma2_SDll_r sigma2_pred_d sigma2_ll95_d sigma2_ul95_d, j(time) i(count)"
global graphhome "Bootstrap-Simulation/fall-2023-RR-submission/figures"


* -----------------
* Create program to re-load sims for each scenario to avoid using too many preserve/restore commands:
cap prog drop reload
prog define reload 
use "data/hellwig/ukdata1.dta", clear
keep timevar fdsupp84 fdunem33 fdinfl21 fdpexp17 polevent falkland el79 el83 el87 el92 el97 trade wilscall thatcher blair // keep only used vars
compress
set seed 03858420

* merge in all the betas across the 3 approaches:
set obs 1001
gen bootcount = _n // to merge on
joinby bootcount using "data/hellwig/betas-parametric.dta"
forv i = 1/22 {
	rename b`i' b_p`i'
}
joinby bootcount using "data/hellwig/betas-me.dta"
forv i = 1/22 {
	rename b`i' b_m`i'
}
joinby bootcount using "data/hellwig/betas-resid.dta"
forv i = 1/22 {
	rename b`i' b_r`i'
}

* b_p = parametric, b_m = maximum entropy, b_r = residual
* We'll create 6 measures of uncertainty: percentiles and +/- SD*1.96 (2 each across the 4 approaches). The line will be the estimated sigma^2 from the OG model for SD method, and median for percentiles

* Run the OG model and save estimates:
arch fdsupp84 fdunem33 fdinfl21 fdpexp17 polevent falkland el79 el83 el87 el92 el97 trade wilscall thatcher blair, robust nolog arch(1) garch(1) het(trade wilscall thatcher blair)
est sto garchmodel
mat B = e(b)
mat list B
* save true B's for type 2 CI creation
svmat B, names(trueB_)
forv i = 1/22 {
	su trueB_`i' in 1
	replace trueB_`i' = r(mean)
}
/* here's the important ones, the HET equation:
beta_16 = trade
beta_17 = wilscall
beta_18 = thatcher
beta_19 = blair 
beta_20 = _cons 
beta_21 = arch
beta_22 = garch
*/
end
* -----------------




* ----------------------------------------------------------------------------
*	Table 2 (main manuscript): Trade exposure and UK Government Support, Hellwig (2007)
* ----------------------------------------------------------------------------
reload // GARCH regression table replicates Table 2
* ----------------------------------------------------------------------------






reload // re-load sims for this scenario




* ----------------------------------------------------------------------------
*	Trade shock (+4points) under WILSON/CALLAHAN
*	[This code replicates all the + trade shock Wilson/Callahan
*	 plots in the manuscript and SI]
* ----------------------------------------------------------------------------
* wilson/callahan is beta_17
* First create OG predictions:
su trade, meanonly
loc tradeval = r(mean)
gen sigma2true_0 = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*1 + trueB_18*0 + trueB_19*0) // this is the one that we're going to use as our main prediction 
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*1 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)
}
* shock: at t=31, shock tradeval by mean + 4:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_20 + trueB_16*(`tradeval'+ 4) + trueB_17*1 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)

* post-shock, at t = 32+, put tradeval to mean and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*(`tradeval') + trueB_17*1 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
est replay garchmodel
su trade, meanonly
loc tradeval = r(mean)
* init
nlcom exp(_b[HET:_cons] + _b[HET:trade]*`tradeval' + _b[HET:wilscall]*1 + _b[HET:thatcher]*0 + _b[HET:blair]*0) 
mat delta_b = r(table)
* now loop t=1/30:
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*1 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t=31, shock tradeval by mean + 4:
nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*(`tradeval'+4) + _b[HET:wilscall]*1 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
mat delta_b = r(table)
gen sigma2_pred_d31 = delta_b[1,1] // expected value
gen sigma2_ll95_d31 = delta_b[5,1] // ll
gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, tradeval back at mean and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*1 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r
	su trade, meanonly
	loc tradeval = r(mean)
	gen sigma2_`b'0 = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*1 + b_`b'18*0 + b_`b'19*0)
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*1 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)
	}
	* shock: at t=31, shock tradeval by mean + 4:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'20 + b_`b'16*(`tradeval'+ 4) + b_`b'17*1 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)

	* post-shock, at t = 32+, keep tradeval at mean and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*(`tradeval') + b_`b'17*1 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
keep in 29/40
drop time
gen time = _n - 1
	
	* Parametric percentile (part of Wilson/Callahan row in FIgure SM.11)
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/tradewillcall-parametric-percentile.pdf", as(pdf) replace
	* ME percentile (part of Wilson/Callahan row in FIgure SM.11)
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/tradewillcall-me-percentile.pdf", as(pdf) replace
	* resid percentile (part of Wilson/Callahan row in FIgure SM.11)
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/tradewillcall-resid-percentile.pdf", as(pdf) replace
	
	
	* now parametric SD (part of Wilson/Callahan row in Figure SM.12)
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/tradewillcall-parametric-sd.pdf", as(pdf) replace
	* ME SD (part of Wilson/Callahan row in Figure SM.12)
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/tradewillcall-me-sd.pdf", as(pdf) replace
	* resid SD (part of Wilson/Callahan row in Figure SM.12)
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/tradewillcall-resid-sd.pdf", as(pdf) replace
	
	* delta (Figure SM.1 (a))
	twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs and expected values from delta method") title("Wilson/Callahan")
	*graph export "$graphhome/tradewillcall-delta-sd.pdf", as(pdf) replace

* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}	
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

	* Parametric percentile (part of Wilson/Callahan row in Figure 2)
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Wilson/Callahan") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/tradewillcall-parametric-percentile-rescale.pdf", as(pdf) replace
	* ME percentile (part of Wilson/Callahan row in Figure 2)
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Wilson/Callahan") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/tradewillcall-me-percentile-rescale.pdf", as(pdf) replace
	* resid percentile (part of Wilson/Callahan row in Figure 2)
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Wilson/Callahan") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/tradewillcall-resid-percentile-rescale.pdf", as(pdf) replace
	
	
	* delta: (Figure 1 (b))
	twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method") title("Wilson/Callahan")
	*graph export "$graphhome/tradewillcall-delta-percentile-rescale.pdf", as(pdf) replace


* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad

* Delta (part of the Wilson/Callahan row in Figure SM.13)
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method") title("Wilson/Callahan")
	*graph export "$graphhome/tradewillcall-delta-effectsplot.pdf", as(pdf) replace

* parametric (part of the Wilson/Callahan row in Figure SM.13)
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles") title("Wilson/Callahan")
	*graph export "$graphhome/tradewillcall-parametric-effectsplot.pdf", as(pdf) replace
	
* ME (part of the Wilson/Callahan row in Figure SM.13)
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles") title("Wilson/Callahan")
	*graph export "$graphhome/tradewillcall-me-effectsplot.pdf", as(pdf) replace
	
* resid (part of the Wilson/Callahan row in Figure SM.13)
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles") title("Wilson/Callahan")
	*graph export "$graphhome/tradewillcall-resid-effectsplot.pdf", as(pdf) replace	
* ----------------------------------------------------------------------------




reload // re-load sims for this scenario





* ----------------------------------------------------------------------------
*	Trade shock (+4points) under THATCHER
*	[This code replicates all the + trade shock Thatcher
*	 plots in the manuscript and SI]
* ----------------------------------------------------------------------------
* thathcer is beta_18
* First create OG predictions:
su trade, meanonly
loc tradeval = r(mean)
gen sigma2true_0 = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*1 + trueB_19*0) // this is the one that we're going to use as our main prediction 
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*1 + trueB_19*0) + trueB_22*r(mean)
}
* shock: at t=31, shock tradeval by mean + 4:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_20 + trueB_16*(`tradeval'+ 4) + trueB_17*0 + trueB_18*1 + trueB_19*0) + trueB_22*r(mean)

* post-shock, at t = 32+, keep tradeval at mean and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*(`tradeval') + trueB_17*0 + trueB_18*1 + trueB_19*0) + trueB_22*r(mean)
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
su trade, meanonly
loc tradeval = r(mean)
* init
nlcom exp(_b[HET:_cons] + _b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + _b[HET:thatcher]*1 + _b[HET:blair]*0) 
mat delta_b = r(table)
* now loop t=1/30:
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*1 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t=31, shock tradeval by mean + 4:
nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*(`tradeval'+4) + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*1 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
mat delta_b = r(table)
gen sigma2_pred_d31 = delta_b[1,1] // expected value
gen sigma2_ll95_d31 = delta_b[5,1] // ll
gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, tradeval back at mean and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*1 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r
	su trade, meanonly
	loc tradeval = r(mean)
	gen sigma2_`b'0 = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*1 + b_`b'19*0)
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*1 + b_`b'19*0) + b_`b'22*r(mean)
	}
	* shock: at t=31, shock tradeval by mean + 4:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'20 + b_`b'16*(`tradeval'+ 4) + b_`b'17*0 + b_`b'18*1 + b_`b'19*0) + b_`b'22*r(mean)

	* post-shock, at t = 32+, keep tradeval at mean and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*(`tradeval') + b_`b'17*0 + b_`b'18*1 + b_`b'19*0) + b_`b'22*r(mean)
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
keep in 29/40
drop time
gen time = _n - 1
	
* each plot separate
	* Parametric percentile (part of Thatcher row in Figure SM.11)
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/tradethatcher-parametric-percentile.pdf", as(pdf) replace
	* ME percentile (part of Thatcher row in Figure SM.11)
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/tradethatcher-me-percentile.pdf", as(pdf) replace
	* resid percentile (part of Thatcher row in Figure SM.11)
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are residual bootstrap percentiles")

	
	* now parametric SD (part of Thatcher row in Figure SM.12)
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/tradethatcher-parametric-sd.pdf", as(pdf) replace
	* ME SD (part of Thatcher row in Figure SM.12)
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/tradethatcher-me-sd.pdf", as(pdf) replace
	* resid SD (part of Thatcher row in Figure SM.12)
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/tradethatcher-resid-sd.pdf", as(pdf) replace
	
	* delta (Figure SM.1 (c))
	twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs and expected values from delta method")
	*graph export "$graphhome/tradethatcher-delta-sd.pdf", as(pdf) replace

	
* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

* Parametric percentile (part of Thatcher row in Figure 2)
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Thatcher") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/tradethatcher-parametric-percentile-rescale.pdf", as(pdf) replace
	* ME percentile (part of Thatcher row in Figure 2)
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Thatcher") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/tradethatcher-me-percentile-rescale.pdf", as(pdf) replace
	* resid percentile (part of Thatcher row in Figure 2)
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Thatcher") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/tradethatcher-resid-percentile-rescale.pdf", as(pdf) replace
	

	* do delta: (Figure SM.1 (d))
	twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method")  title("Thatcher")
	*graph export "$graphhome/tradethatcher-delta-percentile-rescale.pdf", as(pdf) replace
	

* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad

* Delta (part of Thatcher row in Figure SM.13)
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method") title("Thatcher")
	*graph export "$graphhome/tradethatcher-delta-effectsplot.pdf", as(pdf) replace

* parametric  (part of Thatcher row in Figure SM.13)
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles") title("Thatcher")
	*graph export "$graphhome/tradethatcher-parametric-effectsplot.pdf", as(pdf) replace
	
* ME  (part of Thatcher row in Figure SM.13)
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles") title("Thatcher")
	*graph export "$graphhome/tradethatcher-me-effectsplot.pdf", as(pdf) replace
	
* resid  (part of Thatcher row in Figure SM.13)
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles") title("Thatcher")
	*graph export "$graphhome/tradethatcher-resid-effectsplot.pdf", as(pdf) replace
* ----------------------------------------------------------------------------




reload // re-load sims for this scenario




* ----------------------------------------------------------------------------
*	Trade shock (+4points) under MAJOR
*	[This code replicates all the + trade shock under Major]
* ----------------------------------------------------------------------------
* major is omitted category
* First create OG predictions:
su trade, meanonly
loc tradeval = r(mean)
gen sigma2true_0 = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*0 + trueB_19*0) // this is the one that we're going to use as our main prediction 
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)
}
* shock: at t=31, shock tradeval by mean + 4:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_20 + trueB_16*(`tradeval'+ 4) + trueB_17*0 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)

* post-shock, at t = 32+, keep tradeval at mean and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*(`tradeval') + trueB_17*0 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
su trade, meanonly
loc tradeval = r(mean)
* init
nlcom exp(_b[HET:_cons] + _b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + _b[HET:thatcher]*0 + _b[HET:blair]*0) 
mat delta_b = r(table)
* now loop t=1/30:
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t=31, shock tradeval by mean + 4:
nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*(`tradeval'+4) + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
mat delta_b = r(table)
gen sigma2_pred_d31 = delta_b[1,1] // expected value
gen sigma2_ll95_d31 = delta_b[5,1] // ll
gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, tradeval back at mean and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r
	su trade, meanonly
	loc tradeval = r(mean)
	gen sigma2_`b'0 = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*0 + b_`b'19*0)
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)
	}
	* shock: at t=31, shock tradeval by mean + 4:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'20 + b_`b'16*(`tradeval'+ 4) + b_`b'17*0 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)

	* post-shock, at t = 32+, keep tradeval at mean and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*(`tradeval') + b_`b'17*0 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
keep in 29/40
drop time
gen time = _n - 1
	
* each plot separate
	* Parametric percentile (part of Major row in Figure SM.11)
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/trademajor-parametric-percentile.pdf", as(pdf) replace
	* ME percentile (part of Major row in Figure SM.11)
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/trademajor-me-percentile.pdf", as(pdf) replace
	* resid percentile (part of Major row in Figure SM.11)
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/trademajor-resid-percentile.pdf", as(pdf) replace

	
	* now parametric SD (part of Major row in Figure SM.12)
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/trademajor-parametric-sd.pdf", as(pdf) replace
	* ME SD (part of Major row in Figure SM.12)
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/trademajor-me-sd.pdf", as(pdf) replace
	* resid SD (part of Major row in Figure SM.12)
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/trademajor-resid-sd.pdf", as(pdf) replace
	
	* delta (Figure SM.1 (e))
	twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs and expected values from delta method")
	*graph export "$graphhome/trademajor-delta-sd.pdf", as(pdf) replace
	
	
* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

* Parametric percentile (part of Major row in Figure 2)
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance""(as % of Pre-Shock Variance)") title("Major") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/trademajor-parametric-percentile-rescale.pdf", as(pdf) replace
	* ME percentile (part of Major row in Figure 2)
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Major") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/trademajor-me-percentile-rescale.pdf", as(pdf) replace
	* resid percentile (part of Major row in Figure 2)
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Major") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/trademajor-resid-percentile-rescale.pdf", as(pdf) replace
	
	* do delta: (Figure SM.1 (f))
	twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method") title("Major") 
	*graph export "$graphhome/trademajor-delta-percentile-rescale.pdf", as(pdf) replace
	
* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad

* Delta (part of Major row in Figure SM.13)
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method") title("Major")
	*graph export "$graphhome/trademajor-delta-effectsplot.pdf", as(pdf) replace

* parametric (part of Major row in Figure SM.13)
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles") title("Major")
	*graph export "$graphhome/trademajor-parametric-effectsplot.pdf", as(pdf) replace
	
* ME (part of Major row in Figure SM.13)
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles") title("Major")
	*graph export "$graphhome/trademajor-me-effectsplot.pdf", as(pdf) replace
	
* resid (part of Major row in Figure SM.13)
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles") title("Major")
	*graph export "$graphhome/trademajor-resid-effectsplot.pdf", as(pdf) replace
* ----------------------------------------------------------------------------







reload // re-load sims for this scenario

* ----------------------------------------------------------------------------
*	Trade shock (+4points) under BLAIR 
*	[This code replicates all the + trade shock Blair plots in the manuscript and SI]
* ----------------------------------------------------------------------------
* First create OG predictions
su trade, meanonly
loc tradeval = r(mean)
gen sigma2true_0 = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*0 + trueB_19*1) // this is the one that we're going to use as our main prediction 
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*0 + trueB_19*1) + trueB_22*r(mean)
}
* shock: at t=31, shock tradeval by mean + 4:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_20 + trueB_16*(`tradeval'+ 4) + trueB_17*0 + trueB_18*0 + trueB_19*1) + trueB_22*r(mean)

* post-shock, at t = 32+, keep tradeval at mean and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*(`tradeval') + trueB_17*0 + trueB_18*0 + trueB_19*1) + trueB_22*r(mean)
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
su trade, meanonly
loc tradeval = r(mean)
* init
nlcom exp(_b[HET:_cons] + _b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + _b[HET:thatcher]*0 + _b[HET:blair]*1) 
mat delta_b = r(table)
* now loop t=1/30:
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*1) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t=31, shock tradeval by mean + 4:
nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*(`tradeval'+4) + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*1) + _b[ARCH:L.arch]*0
mat delta_b = r(table)
gen sigma2_pred_d31 = delta_b[1,1] // expected value
gen sigma2_ll95_d31 = delta_b[5,1] // ll
gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, tradeval back at mean and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*1) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

	
* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r 
	su trade, meanonly
	loc tradeval = r(mean)
	gen sigma2_`b'0 = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*0 + b_`b'19*1)
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*0 + b_`b'19*1) + b_`b'22*r(mean)
	}
	* shock: at t=31, shock tradeval by mean + 4:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'20 + b_`b'16*(`tradeval'+ 4) + b_`b'17*0 + b_`b'18*0 + b_`b'19*1) + b_`b'22*r(mean)

	* post-shock, at t = 32+, put tradeval to mean and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*(`tradeval') + b_`b'17*0 + b_`b'18*0 + b_`b'19*1) + b_`b'22*r(mean)
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
keep in 29/40
drop time
gen time = _n - 1
	
	
	
* each plot is separately created then combined in .tex:

	* Parametric percentile (part of Blair row in Figure SM.11)
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/tradeblair-parametric-percentile.pdf", as(pdf) replace
	* ME percentile (part of Blair row in Figure SM.11)
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/tradeblair-me-percentile.pdf", as(pdf) replace
	* resid percentile (part of Blair row in Figure SM.11)
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/tradeblair-resid-percentile.pdf", as(pdf) replace
	
	
	* now parametric SD (part of Blair row in Figure SM.12)
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/tradeblair-parametric-sd.pdf", as(pdf) replace
	* ME SD (part of Blair row in Figure SM.12)
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/tradeblair-me-sd.pdf", as(pdf) replace
	* resid SD (part of Blair row in Figure SM.12)
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/tradeblair-resid-sd.pdf", as(pdf) replace
	
	* Figure SM.1(g)
	* delta (for delta we'll show the underlying delta expectation)
	twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs and expected values from delta method")
	*graph export "$graphhome/tradeblair-delta-sd.pdf", as(pdf) replace

	
* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}	
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

	* Parametric percentile (part of Blair row in Figure 2)
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Blair") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/tradeblair-parametric-percentile-rescale.pdf", as(pdf) replace
	* ME percentile (part of Blair row in Figure 2)
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Blair") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/tradeblair-me-percentile-rescale.pdf", as(pdf) replace
	* resid percentile (part of Blair row in Figure 2)
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Blair") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/tradeblair-resid-percentile-rescale.pdf", as(pdf) replace
	
	
	* delta: (Figure SM.1 (h))
	twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method") title("Blair")
	*graph export "$graphhome/tradeblair-delta-percentile-rescale.pdf", as(pdf) replace
	
	
* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad

* Delta (part of Blair row in Figure SM.13)
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method") title("Blair")
	*graph export "$graphhome/tradeblair-delta-effectsplot.pdf", as(pdf) replace

* parametric (part of Blair row in Figure SM.13)
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles") title("Blair")
	*graph export "$graphhome/tradeblair-parametric-effectsplot.pdf", as(pdf) replace
	
* ME (part of Blair row in Figure SM.13)
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles") title("Blair")
	*graph export "$graphhome/tradeblair-me-effectsplot.pdf", as(pdf) replace
	
* resid (part of Blair row in Figure SM.13)
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles") title("Blair")
	*graph export "$graphhome/tradeblair-resid-effectsplot.pdf", as(pdf) replace

* ----------------------------------------------------------------------------








reload // re-load sims for this scenario










* ----------------------------------------------------------------------------
*	Trade shock (-4points) under WILSON/CALLAHAN
*	[This code replicates all the - trade shock Wilson/Callahan
*	 plots in the manuscript and SI]
* ----------------------------------------------------------------------------
* wilson/callahan is beta_17
* First create OG predictions:
su trade, meanonly
loc tradeval = r(mean)
gen sigma2true_0 = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*1 + trueB_18*0 + trueB_19*0) // this is the one that we're going to use as our main prediction 
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*1 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)
}
* shock: at t=31, shock tradeval by mean - 4:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_20 + trueB_16*(`tradeval'- 4) + trueB_17*1 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)

* post-shock, at t = 32+, keep tradeval at mean and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*(`tradeval') + trueB_17*1 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
su trade, meanonly
loc tradeval = r(mean)
* init
nlcom exp(_b[HET:_cons] + _b[HET:trade]*`tradeval' + _b[HET:wilscall]*1 + _b[HET:thatcher]*0 + _b[HET:blair]*0) 
mat delta_b = r(table)
* now loop t=1/30:
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*1 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t=31, shock tradeval by mean + 4:
nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*(`tradeval'-4) + _b[HET:wilscall]*1 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
mat delta_b = r(table)
gen sigma2_pred_d31 = delta_b[1,1] // expected value
gen sigma2_ll95_d31 = delta_b[5,1] // ll
gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, tradeval back at mean and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*1 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r 
	su trade, meanonly
	loc tradeval = r(mean)
	gen sigma2_`b'0 = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*1 + b_`b'18*0 + b_`b'19*0)
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*1 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)
	}
	* shock: at t=31, shock tradeval by mean - 4:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'20 + b_`b'16*(`tradeval'- 4) + b_`b'17*1 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)

	* post-shock, at t = 32+, keep tradeval at mean and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*(`tradeval') + b_`b'17*1 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
keep in 29/40
drop time
gen time = _n - 1
	
* each plot separate
	* Parametric percentile (part of Wilson/Callahan row in Figure SM.14)
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/negtradewillcall-parametric-percentile.pdf", as(pdf) replace
	* ME percentile (part of Wilson/Callahan row in Figure SM.14)
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/negtradewillcall-me-percentile.pdf", as(pdf) replace
	* resid percentile (part of Wilson/Callahan row in Figure SM.14)
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/negtradewillcall-resid-percentile.pdf", as(pdf) replace

	
	* now parametric SD (part of Wilson/Callahan row in Figure SM.15)
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/negtradewillcall-parametric-sd.pdf", as(pdf) replace
	* ME SD  (part of Wilson/Callahan row in Figure SM.15)
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/negtradewillcall-me-sd.pdf", as(pdf) replace
	* resid SD  (part of Wilson/Callahan row in Figure SM.15)
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/negtradewillcall-resid-sd.pdf", as(pdf) replace
	
	* delta (Figure SM.2 (a))
	twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Wilson/Callahan") note("CIs and expected values from delta method")
	*graph export "$graphhome/negtradewillcall-delta-sd.pdf", as(pdf) replace


* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

	* Parametric percentile (part of Wilson/Callahan row in Figure 3)
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Wilson/Callahan") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/negtradewillcall-parametric-percentile-rescale.pdf", as(pdf) replace
	* ME percentile (part of Wilson/Callahan row in Figure 3)
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Wilson/Callahan") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/negtradewillcall-me-percentile-rescale.pdf", as(pdf) replace
	* resid percentile (part of Wilson/Callahan row in Figure 3)
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Wilson/Callahan") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/negtradewillcall-resid-percentile-rescale.pdf", as(pdf) replace

	
	* delta: (Figure SM.2 (b))
	twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method") title("Wilson/Callahan")
	*graph export "$graphhome/negtradewillcall-delta-percentile-rescale.pdf", as(pdf) replace
	
* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad

* parametric (part of Wilson/Callahan row in Figure SM.16)
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles") title("Wilson/Callahan")
	*graph export "$graphhome/negtradewillcall-parametric-effectsplot.pdf", as(pdf) replace
	
* ME (part of Wilson/Callahan row in Figure SM.16)
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles") title("Wilson/Callahan")
	*graph export "$graphhome/negtradewillcall-me-effectsplot.pdf", as(pdf) replace
	
* resid (part of Wilson/Callahan row in Figure SM.16)
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles") title("Wilson/Callahan")
	*graph export "$graphhome/negtradewillcall-resid-effectsplot.pdf", as(pdf) replace

* Delta (part of Wilson/Callahan row in Figure SM.16)
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method") title("Wilson/Callahan")
	*graph export "$graphhome/negtradewillcall-delta-effectsplot.pdf", as(pdf) replace
* ----------------------------------------------------------------------------





reload // re-load sims for this scenario




* ----------------------------------------------------------------------------
*	Trade shock (-4points) under THATCHER
*	[THis code replicates all the - trade shock Thatcher plots in
*	the manuscript and SI]
* ----------------------------------------------------------------------------
* thathcer is beta_18
* First create OG predictions:
su trade, meanonly
loc tradeval = r(mean)
gen sigma2true_0 = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*1 + trueB_19*0) // this is the one that we're going to use as our main prediction 
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*1 + trueB_19*0) + trueB_22*r(mean)
}
* shock: at t=31, shock tradeval by mean - 4:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_20 + trueB_16*(`tradeval'- 4) + trueB_17*0 + trueB_18*1 + trueB_19*0) + trueB_22*r(mean)

* post-shock, at t = 32+, keep tradeval at mean and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*(`tradeval') + trueB_17*0 + trueB_18*1 + trueB_19*0) + trueB_22*r(mean)
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
su trade, meanonly
loc tradeval = r(mean)
* init
nlcom exp(_b[HET:_cons] + _b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + _b[HET:thatcher]*1 + _b[HET:blair]*0) 
mat delta_b = r(table)
* now loop t=1/30:
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*1 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t=31, shock tradeval by mean + 4:
nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*(`tradeval'-4) + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*1 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
mat delta_b = r(table)
gen sigma2_pred_d31 = delta_b[1,1] // expected value
gen sigma2_ll95_d31 = delta_b[5,1] // ll
gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, tradeval back at mean and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*1 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r
	su trade, meanonly
	loc tradeval = r(mean)
	gen sigma2_`b'0 = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*1 + b_`b'19*0)
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*1 + b_`b'19*0) + b_`b'22*r(mean)
	}
	* shock: at t=31, shock tradeval by mean - 4:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'20 + b_`b'16*(`tradeval'- 4) + b_`b'17*0 + b_`b'18*1 + b_`b'19*0) + b_`b'22*r(mean)

	* post-shock, at t = 32+, keep tradeval at mean and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*(`tradeval') + b_`b'17*0 + b_`b'18*1 + b_`b'19*0) + b_`b'22*r(mean)
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
keep in 29/40
drop time
gen time = _n - 1
	
* each plot separate
	* Parametric percentile (part of Thatcher row in Figure SM.14)
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/negtradethatcher-parametric-percentile.pdf", as(pdf) replace
	* ME percentile (part of Thatcher row in Figure SM.14)
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/negtradethatcher-me-percentile.pdf", as(pdf) replace
	* resid percentile (part of Thatcher row in Figure SM.14)
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/negtradethatcher-resid-percentile.pdf", as(pdf) replace
	
	
	* now parametric SD (part of Thatcher row in Figure SM.15)
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/negtradethatcher-parametric-sd.pdf", as(pdf) replace
	* ME SD (part of Thatcher row in Figure SM.15)
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/negtradethatcher-me-sd.pdf", as(pdf) replace
	* resid SD (part of Thatcher row in Figure SM.15)
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/negtradethatcher-resid-sd.pdf", as(pdf) replace
	
	* delta (Figure SM.2 (c))
	twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Thatcher") note("CIs and expected values from delta method")
	*graph export "$graphhome/negtradethatcher-delta-sd.pdf", as(pdf) replace

	
* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

	* Parametric percentile (part of Thatcher row in Figure 3)
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Thatcher") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/negtradethatcher-parametric-percentile-rescale.pdf", as(pdf) replace
	* ME percentile (part of Thatcher row in Figure 3)
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Thatcher") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/negtradethatcher-me-percentile-rescale.pdf", as(pdf) replace
	* resid percentile (part of Thatcher row in Figure 3)
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Thatcher") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/negtradethatcher-resid-percentile-rescale.pdf", as(pdf) replace

	
	* delta: (Figure SM.2 (d))
	twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method") title("Thatcher")
	*graph export "$graphhome/negtradethatcher-delta-percentile-rescale.pdf", as(pdf) replace
	
* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad

* parametric  (part of Thatcher row in Figure SM.16)
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles") title("Thatcher")
	*graph export "$graphhome/negtradethatcher-parametric-effectsplot.pdf", as(pdf) replace
	
* ME  (part of Thatcher row in Figure SM.16)
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles") title("Thatcher")
	*graph export "$graphhome/negtradethatcher-me-effectsplot.pdf", as(pdf) replace
	
* resid  (part of Thatcher row in Figure SM.16)
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles") title("Thatcher")
	*graph export "$graphhome/negtradethatcher-resid-effectsplot.pdf", as(pdf) replace

	* Delta (part of Thatcher row in Figure SM.16)
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method") title("Thatcher")
	*graph export "$graphhome/negtradethatcher-delta-effectsplot.pdf", as(pdf) replace
* ----------------------------------------------------------------------------





reload // re-load sims for this scenario



* ----------------------------------------------------------------------------
*	Trade shock (-4points) under MAJOR
*	[this code replicates all the - trade shock Major plots in the manuscript and SI]
* ----------------------------------------------------------------------------
* major is omitted category
* First create OG predictions:
su trade, meanonly
loc tradeval = r(mean)
gen sigma2true_0 = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*0 + trueB_19*0) // this is the one that we're going to use as our main prediction 
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)
}
* shock: at t=31, shock tradeval by mean - 4:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_20 + trueB_16*(`tradeval'- 4) + trueB_17*0 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)

* post-shock, at t = 32+, keep tradeval at mean and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*(`tradeval') + trueB_17*0 + trueB_18*0 + trueB_19*0) + trueB_22*r(mean)
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
su trade, meanonly
loc tradeval = r(mean)
* init
nlcom exp(_b[HET:_cons] + _b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + _b[HET:thatcher]*0 + _b[HET:blair]*0) 
mat delta_b = r(table)
* now loop t=1/30:
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t=31, shock tradeval by mean + 4:
nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*(`tradeval'-4) + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
mat delta_b = r(table)
gen sigma2_pred_d31 = delta_b[1,1] // expected value
gen sigma2_ll95_d31 = delta_b[5,1] // ll
gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, tradeval back at mean and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*0) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r 
	su trade, meanonly
	loc tradeval = r(mean)
	gen sigma2_`b'0 = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*0 + b_`b'19*0)
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)
	}
	* shock: at t=31, shock tradeval by mean - 4:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'20 + b_`b'16*(`tradeval'- 4) + b_`b'17*0 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)

	* post-shock, at t = 32+, keep tradeval at mean and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*(`tradeval') + b_`b'17*0 + b_`b'18*0 + b_`b'19*0) + b_`b'22*r(mean)
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
keep in 29/40
drop time
gen time = _n - 1
	
* each plot separate
	* Parametric percentile  (part of Major row in Figure SM.14)
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/negtrademajor-parametric-percentile.pdf", as(pdf) replace
	* ME percentile (part of Major row in Figure SM.14)
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/negtrademajor-me-percentile.pdf", as(pdf) replace
	* resid percentile (part of Major row in Figure SM.14)
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/negtrademajor-resid-percentile.pdf", as(pdf) replace
	
	
	* now parametric SD (part of Major row in Figure SM.15)
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/negtrademajor-parametric-sd.pdf", as(pdf) replace
	* ME SD  (part of Major row in Figure SM.15)
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/negtrademajor-me-sd.pdf", as(pdf) replace
	* resid SD  (part of Major row in Figure SM.15)
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/negtrademajor-resid-sd.pdf", as(pdf) replace
	
	* delta (Figure SM.2 (e))
	twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Major") note("CIs and expected values from delta method")
	*graph export "$graphhome/negtrademajor-delta-sd.pdf", as(pdf) replace
	

* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

	* Parametric percentile  (part of Major row in Figure 3)
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Major") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/negtrademajor-parametric-percentile-rescale.pdf", as(pdf) replace
	* ME percentile (part of Major row in Figure 3)
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Major") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/negtrademajor-me-percentile-rescale.pdf", as(pdf) replace
	* resid percentile (part of Major row in Figure 3)
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Major") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/negtrademajor-resid-percentile-rescale.pdf", as(pdf) replace
	
	* delta: (Figure SM.2 (h))
	twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method") title("Major")
	*graph export "$graphhome/negtrademajor-delta-percentile-rescale.pdf", as(pdf) replace
	
* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad

* parametric (part of Major row in Figure SM.16)
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles") title("Major")
	*graph export "$graphhome/negtrademajor-parametric-effectsplot.pdf", as(pdf) replace
	
* ME (part of Major row in Figure SM.16)
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles") title("Major")
	*graph export "$graphhome/negtrademajor-me-effectsplot.pdf", as(pdf) replace
	
* resid (part of Major row in Figure SM.16)
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles") title("Major")
	*graph export "$graphhome/negtrademajor-resid-effectsplot.pdf", as(pdf) replace

* Delta (part of Major row in Figure SM.16)
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method") title("Major")
	*graph export "$graphhome/negtrademajor-delta-effectsplot.pdf", as(pdf) replace
* ----------------------------------------------------------------------------







reload // re-load sims for this scenario




* ----------------------------------------------------------------------------
*	Trade shock (-4points) under BLAIR
*	[This code replicates all the - trade shock Blair plots in the manuscript and SI]
* ----------------------------------------------------------------------------
* First create OG predictions:
su trade, meanonly
loc tradeval = r(mean)
gen sigma2true_0 = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*0 + trueB_19*1) // this is the one that we're going to use as our main prediction 
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*`tradeval' + trueB_17*0 + trueB_18*0 + trueB_19*1) + trueB_22*r(mean)
}
* shock: at t=31, shock tradeval by mean + 4:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_20 + trueB_16*(`tradeval'- 4) + trueB_17*0 + trueB_18*0 + trueB_19*1) + trueB_22*r(mean)

* post-shock, at t = 32+, keep tradeval at mean and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_20 + trueB_16*(`tradeval') + trueB_17*0 + trueB_18*0 + trueB_19*1) + trueB_22*r(mean)
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
su trade, meanonly
loc tradeval = r(mean)
* init
nlcom exp(_b[HET:_cons] + _b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + _b[HET:thatcher]*0 + _b[HET:blair]*1) 
mat delta_b = r(table)
* now loop t=1/30:
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*1) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t=31, shock tradeval by mean + 4:
nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*(`tradeval'-4) + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*1) + _b[ARCH:L.arch]*0
mat delta_b = r(table)
gen sigma2_pred_d31 = delta_b[1,1] // expected value
gen sigma2_ll95_d31 = delta_b[5,1] // ll
gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, tradeval back at mean and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + exp(_b[HET:_cons] + ///
		_b[HET:trade]*`tradeval' + _b[HET:wilscall]*0 + ///
		_b[HET:thatcher]*0 + _b[HET:blair]*1) + _b[ARCH:L.arch]*0
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r 
	su trade, meanonly
	loc tradeval = r(mean)
	gen sigma2_`b'0 = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*0 + b_`b'19*1)
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*`tradeval' + b_`b'17*0 + b_`b'18*0 + b_`b'19*1) + b_`b'22*r(mean)
	}
	* shock: at t=31, shock tradeval by mean + 4:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'20 + b_`b'16*(`tradeval'- 4) + b_`b'17*0 + b_`b'18*0 + b_`b'19*1) + b_`b'22*r(mean)

	* post-shock, at t = 32+, keep tradeval at mean and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'20 + b_`b'16*(`tradeval') + b_`b'17*0 + b_`b'18*0 + b_`b'19*1) + b_`b'22*r(mean)
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
keep in 29/40
drop time
gen time = _n - 1
	
* each plot separate
	* Parametric percentile (part of Blair row in Figure SM.14)
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/negtradeblair-parametric-percentile.pdf", as(pdf) replace
	* ME percentile (part of Blair row in Figure SM.14)
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/negtradeblair-me-percentile.pdf", as(pdf) replace
	* resid percentile (part of Blair row in Figure SM.14)
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/negtradeblair-resid-percentile.pdf", as(pdf) replace
	
	* now parametric SD (part of Blair row in Figure SM.15)
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/negtradeblair-parametric-sd.pdf", as(pdf) replace
	* ME SD (part of Blair row in Figure SM.15)
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/negtradeblair-me-sd.pdf", as(pdf) replace
	* resid SD (part of Blair row in Figure SM.15)
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/negtradeblair-resid-sd.pdf", as(pdf) replace
	
	* delta (Figure SM.2 (g))
	twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") title("Blair") note("CIs and expected values from delta method")
	*graph export "$graphhome/negtradeblair-delta-sd.pdf", as(pdf) replace


	* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

	* Parametric percentile (part of Blair row in Figure 3)
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Blair") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/negtradeblair-parametric-percentile-rescale.pdf", as(pdf) replace
	* ME percentile (part of Blair row in Figure 3)
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Blair") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/negtradeblair-me-percentile-rescale.pdf", as(pdf) replace
	* resid percentile (part of Blair row in Figure 3)
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") title("Blair") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/negtradeblair-resid-percentile-rescale.pdf", as(pdf) replace
	
	* do delta: (Figure SM.2 (h))
	twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method") title("Blair")
	*graph export "$graphhome/negtradeblair-delta-percentile-rescale.pdf", as(pdf) replace
	
* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad

* parametric (part of Blair row in Figure SM.16)
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles") title("Blair")
	*graph export "$graphhome/negtradeblair-parametric-effectsplot.pdf", as(pdf) replace
	
* ME (part of Blair row in Figure SM.16)
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles") title("Blair")
	*graph export "$graphhome/negtradeblair-me-effectsplot.pdf", as(pdf) replace
	
* resid (part of Blair row in Figure SM.16)
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles") title("Blair")
	*graph export "$graphhome/negtradeblair-resid-effectsplot.pdf", as(pdf) replace

* Delta (part of Blair row in Figure SM.16)
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method") title("Blair")
	*graph export "$graphhome/negtradeblair-delta-effectsplot.pdf", as(pdf) replace
* ----------------------------------------------------------------------------




